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We describe standing wave patterns induced by an attractive finite-ranged external potential inside 
a large Bose-Einstein Condensate (BEC). As the potential depth increases, the time independent 
Gross-Pitaevskii equation develops pairs of solutions that have nodes in their wavefunction. We 
elucidate the nature of these states and study their dynamical stability. Although we study the 
problem in a two-dimensional BEC subject to a cylindrically symmetric square well potential of a 
radius that is comparable to the coherence length of the BEC, our analysis reveals general trends, 
valid in two and three dimensions, independent of the symmetry of the localized potential well, 
and suggestive of the behavior in general, short- and large-range potentials. One set of nodal BEC 
wavefunctions resembles the single particle n node bound state wavefunction of the potential well, 
the other wavefunctions resemble the n—1 node bound-state wavefunction with a kink state pinned 
by the potential. The second state, though corresponding to the lower free energy value of the pair 
of n node BEC states, is always unstable, whereas the first can be dynamically stable in intervals 
of the potential well depth, implying that the standing wave BEC can evolve from a dynamically 
unstable to stable, and back to unstable status as the potential well is adiabatically deepened, a 
phenomenon that we refer to as "reentrant dynamical stability" . 

PACS numbers: 03.67.Lx, 73.43. Nq,03. 75. Hh,67.40.Yv 



I. INTRODUCTION 



As a superfluid, the dilute gas Bose-Einstein conden- 
sate (BEC) is a coherent quantum system and behaves in 
many respects like a single-particle quantum wave. Sta- 
tionary single particle bound states have wave functions 
with nodes, giving standing wave patterns. In a BEC, 
nodal structures, such as the one pictured in Fig.[l] could 
be imaged and manipulated directly since the experimen- 
tal BEC wave functions extend over tens of microns. Can 
BECs support such patterns as long-lived structures? A 
fundamental difference between single particle and BEC 
wave dynamics is the nonlinearity of the BEC evolution 
that stems from the inter-particle interactions. This non- 
linearity leads to stationary standing wave solutions such 
as kink states or dark solitons that have been created in 
BEC experiments P, However, the same nonlinear- 
ity also makes the kink states unstable in two and three 
dimensions [MJl, Q , as has been observed in cold atom 
experiments [6[. The work of this paper was motivated 
by the basic question: when and how can a BEC support 
a long-lived standing wave induced by a local potential 
well in dimensions higher than one? The standing wave 
gives a stationary interference pattern in the full single- 
particle density that vanishes [31 on the surfaces at which 
the stationary BEC wavefunction changes sign. The exis- 
tence of such structures may, as wc discuss below, impact 
fundamental science phenomena and cold atom applica- 
tions. 

We analyze the nodal BEC wavefunctions in a two- 
dimensional (pancake shaped) BEC; although the trends 
uncovered by our analysis, suitably generalized, are valid 
in higher dimensions. The 2D trap geometry, realized 
in cold atom experiments Q, offers unprecedented 
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FIG. 1: Example of an atomic nodal BEC density (well radius 
b = 1.0^, well depth 77 = 1.8). This case happens to be in a 
reentrant stable regime. 



prospects for "wavefunction engineering" : focused laser- 
beams can, in principle, access and image each point of 
the BEC wavefunction. We consider a large, dilute 2D 
BEC in a homogeneous potential that is subjected to a 
superimposed attractive potential well of cylindrical sym- 
metry with a radius comparable to the coherence length 
of the surrounding BEC. We study the existence of cylin- 
drically symmetric standing wave solutions of the time- 
independent Gross-Pitaevskii equation, elucidating their 
nature and investigating their dynamical stability. 

Our analysis reveals general trends: as the depth of 
the finite-ranged potential well increases near the reso- 
nance value (though still above) at which a single-particle 
bound state of n nodes forms, a pair of n node standing 
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wave solutions to the time independent Gross-Pitaevskii 
equation appears, as found in Ref. [l0| . Upon further 
deepening of the well, the last node of one of the two 
new BEG wavefunctions moves inside the well and the 
corresponding wavefunction resembles that of the single 
particle bound state. We refer to this solution as the 
"atomic nodal BEG" state. At the same well depth, the 
other n node BEG-wavcfunction resembles the wavefunc- 
tion of the 77,-1 node atomic bound state with an extra 
soliton pinned near the potential edge and we refer to 
it as the "nonlinear nodal BEG" solution. Of the two 
n node BEG wavefunctions, the nonlinear nodal BEG 
state has the lower free energy, but we find that it is 
always dynamically unstable. If the atomic nodal BEG 
wavefunction varies on a length scale comparable to or 
greater than the healing length of the surrounding BEG, 
then the atomic BEG solution is also unstable when it 
first forms. Upon further deepening of the well, the sta- 
tionary atomic nodal BEG wavefunction can become sta- 
ble, then go unstable again as the well depth is further 
increased — a phenomenon that we refer to as "reentrant 
dynamical stability" . We show that the stability analysis 
can be reduced to the simpler study of the spectrum of 
two distinct single-particlc-Hamiltonian-like operators. 

The size of the stability islands, in terms of the range 
of well depths, for the atomic nodal BEG standing wave 
depends sensitively on the ratio of the potential range to 
the BEG coherence length: the more the potential range 
exceeds the coherence length, the broader the regions of 
instability become. If the potential well is confined to 
a spatial region much smaller than the BEG coherence 
length, the regions of instability are confined to narrow 
intervals of the well depth and the time scale at which 
the instability sets in increases markedly. 

The existence of stable, localized BEG standing waves 
could impact many areas of cold atom physics. In cold 
atom traps standing waves could play a role in the physics 
of localized objects moving through a BEG [l3| • If the ob- 
ject interacts strongly and attractively with BEG bosons, 
the dynamics of object acceleration may bring the BEG 
into an unstable nodal state, disturbing the system and 
possibly providing another mechanism of invalidating the 
picture of superfluid, dissipationless motion [isj . A class 
of BEG objects of particular, fundamental interest con- 
sists of a BEG with self-localized polaron-like impurity 
atoms [l^, [l3| • For sufficiently strong impurity-boson at- 
tractions, a standing wave BEG configuration combined 
with a localized, nodal impurity wavefunction may give 
more complex, localized excited state polaron objects. 
Another class of BEG objects that was proposed [l^ 
consists of mcsoscopic ultra-cold molecules. Localized 
potentials of ions embedded in the BEG may capture 
cold atoms in micron-sized orbits. If the BEG can be 
engineered to be in an atomic nodal BEG state, a gentle 
lowering of the boson-boson interactions can bring the 
system into a true molecular state, assuming that no in- 
stabilities are encountered in the adiabatic dynamics. Fi- 
nally, we remark that long-lived, localized standing wave 



BEG patterns may resolve the challenge of observing the 
BEG response to external perturbations. For instance, a 
ring BEG can sense rotation: as the trap rotates in the 
plane of the ring, the BEG cannot follow as such mo- 
tion would imply vorticity. Hence, in the rotating trap 
frame, the BEG fiows. How can we observe that fiow? 
If a local dip in the ring potential supports a standing 
wave BEG with fringes (across which the BEG cannot 
fiow), the rotation would shift the position of the fringes, 
perhaps destroying the standing wave. 

For a more whimsical application of standing wave 
nodal BEG states, one could look to the hypothesized 
dark matter BEGs. If, as has been speculated, the dark 
matter haloes that surround most galaxies are scalar 
BEGs of ultra-light boson particles [ill . [T^ . [Tsj . steep 
gravitational potentials may support standing wave dark 
matter BEG density patterns. 

The relevance of any of the prospects for standing wave 
BEGs hinge on their stability. In this manuscript we 
study the standing wave patterns and their stability in a 
linear stability analysis (dynamical stability). The paper 
is organized as follows: Section II describes the multiple, 
stationary solutions to the Gross-Piteavskii equation in 
the presence of an attractive potential well of finite ra- 
dius. We describe the linear stability analysis in Section 
HI and we analyze the dynamical stability of the standing 
wave solutions in Section IV. In Section V, we describe 
a general analysis based on the study of the spectrum 
of effective single-particle Hamiltonians. Section VI con- 
cludes. 

II. STANDING WAVE BEG WAVEFUNGTIONS 

We consider a stationary, dilute two-dimensional BEG 
of N bosons, each of mass mi,. Such systems can be 
realized with cold atom harmonic oscillator traps by in- 
creasing one of the the trap frequencies lOz to ensure that 
the corresponding energy significantly exceeds the chemi- 
cal potential ^ of the BEG system, huiz ^ The bosons 
interact mutually via a short-range interaction described 
by a contact interaction potential, i;(r — r') ^ gd{r ~ r'), 
where r is a 2D vector and i5(r) the 2D delta function. 
In the dilute regime, the interaction strength g relates 
to the scattering length a that characterizes low energy 
scattering in three dimensions and t o the gr ound state 
extent 1^ of the uj^ frequency, Iz = \Jh|myJJz^ by the re- 
lation g = {\/^h^ /mb){a/lz) [H, [1^. Here we assume 
that the N bosons are contained by a 2D-box potential of 
macroscopic area VL, corresponding to an average density 
TT-o = N /Vl. The BEG is described by a wavefunction (or 
order parameter) V-'o(r) that is normalized by requiring 

d^r|?/;o(r)p = N. In addition, the BEG experiences 
an attractive external potential T4xt(r) of finite range. 
The wavefunction satisfies the time-independent Gross- 
Pitaevskii equation 

(-^V2 + V;xt(r) + <?|^o|')^o = MV'o. (1) 
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Far from the center of the external potential, the wave- 
function tends to a constant value, ip^ — > y^n^, so that the 
chemical potential is related to the density by fi ~ guo- 
The natural length scale of the problem is the coherence 
length of the BEC, ^, defined by h'^/nibS^'^ = fi. Scahng 
the energy by /i, the length by ^, x = r/^, and the density 
by Uo, so that 0o = ipo/ we obtain the dimcnsionlcss 
form of the Gross-Pitaevskii equation, 

(-^ + C/o(x) + |(/)oP)</'o=</'o, (2) 

where Uo{x) denotes the scaled external potential 
f/o(x) = T4xt(r = and the large distance bound- 

ary condition becomes limj;_,oo 0o(x) = 1. 

We perform the calculations for an attractive exter- 
nal potential that has a square well shape of radius b and 
depth TT^fi^T]^ /SrUbb'^. In three dimensions the resonances 
for noninteracting particles in this potential would occur 
when rj is an odd integer; while in 2D the noninteract- 
ing resonances occur at r/ w 2.45, 4.48, . . ., relating to 
solutions of a transcendental equation involving Bessel 
functions. Eq. ([2]) then reads 

(-^- 8^^0(VC--^) + l</'op)</'o = 0o, (3) 

where Q is the unit step function. For a given well radius 
6/^; we vary the dimcnsionlcss well depth j] and solve 
Eq.'®. 

To integrate Eq. Q subject to the proper boundry con- 
ditions, we use a scheme introduced in Ref. [l^l to con- 
struct the time-independent BEC wavefunctions in the 
ion-neutral atom polarization potential. The boundary 
conditions introduce parameters Ai and A2'. at the ori- 
gin, (j)'^{x ^ 0) ^ and (j)o{x ^ 0) = Ai, whereas far 
from the localized potential, (j)o{x) ^ 1 + ^2 exp(— 22;). 
We pick the radius of the potential well, Xm = b/£^, as 
the matching point for the inward and outward integra- 
tion procedure. Integrating inward from large x-values 
for different choices of A2, we find the wavefunction and 
its derivative at the matching point, {(f>o{xm)T 4>'oi^rn)), 
tracing out a parametric curve that is independent of the 
potential (since the equation is integrated in the region 
where L'o(x) = for x > x„i). Similarly, the outward 
integration from a; = for different Ai values yields a 
second {4>oixm),4''oi^m)) curve, this one dependent on 
the potential. As the values of the wavefunction and its 
derivative must match at Xm, the intersections of both 
curves (and the corresponding Ai and A2 parameter val- 
ues) determine the stationary solutions. Fig. [5]shows the 
curves for well radius b = 1.0^ and well depth r/ — 2.9. 
These parameters result in five crossings and thus five 
solutions. As mentioned, the "outside" curve doesn't de- 
pend on well depth rj. When rj is small, there is only 
a single crossing of the outside curve with the "inside" 
curve, corresponding to a single solution — the familiar 
nodeless BEC. As the well depth r] increases, the in- 
side curve starts to swirl and pinches through the out- 
side curve, generating a pair of new solutions. When the 



well depth has reached the 77 = 2.9 value shown in Fig. [2l 
the inside curve pinches once more through the outside 
curve near 0o(&) ~ —0.8, generating two additional time- 
independent Gross-Pitaveskii solutions. 
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FIG. 2: Solving Eq. ((3} graphically, matching inner (solid) 
and outer (dashed) integrations at the potential edge, Xm = 
VC) by locating the crossing points in the {4>'o{b/i,),4>o{h/(,)) 
plane. The example shown, for a well radius of & = 1.0^ and 
well depth 77 = 2.9, gives five crossing points and thus five 
solutions to the GP equation, Eq. (jSj. 

The new solutions differ qualitatively from the nodeless 
BEC solution. As 77 increases, each one of the new pair 
of BEC wavefunctions has one more radial node than 
the number of nodes of the previous pair. At the well 
depth where the new pair of solutions first exists, the 
two new solutions are degenerate and they have the last 
node located slightly beyond the well radius. For a well 
radius of b = 1.0^, wc find the critical well depth for 
the one node solutions is 77 w 1.2444. In Fig. [3l^a) for 
77 — 1.25, the two new solutions have visibly split. As the 
well depth increases, the two nodal solutions continue to 
differentiate: the node of one moves inside the potential 
well as the well depth increases, while the the node of 
the other solution moves outside the potential well as 
depicted in Fig. [3Ub) for well depth ?/ = 2.5. As the 
well depth increases further, two new solutions, each with 
two nodes, appear near 77 sa 2.86; and we show the five 
solutions for 77 = 3.5 in Fig.[3]Jc). 

We clarify the nature of the different nodal solutions. 
In accordance with Levinson's theorem, the deepening 
of the potential well produces bound states in the single 
particle wavefunction description. Each time the zero 
energy scattering phase goes through a resonance (in- 
creasing by 7r), a new bound state forms with one more 
node than the previously formed bound state. We note 
that one class of the time-independent Gross-Pitaevskii 
solutions — the solution in which the last node moves in- 
ward as the potential depth increases — starts to resem- 
ble the noninteracting (NI) zero energy wavefunction. In 
Fig. [3fc) , the solid line curves show three BEC wavefunc- 
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FIG. 3: Wavefunctions ^o/y/n^ = 0o solving Eq. Q for a 
well radius b = 1.0^ (marked as a hash on axis) and well 
depths 77 of (a) 1.25, (b) 2.5, and (c) 3.5. As the well is made 
deeper, multiple solutions are found. 



tions with zero, one, and two nodes (the dashed curve 
wavefunctions will be discussed below). Compare these 
wavefunctions to the NI wavefunction plotted in Fig. 3] 
for the same well depth. The NI bound states with zero 
and one node exhibit a qualitative similarity to the BEC 
wavefunctions with zero and one nodes, except that the 
former vanish at large distance whereas the latter tend 
to a constant value (Figs. [5] and [3Jc)). Moreover, the 
NI zero-energy scattering solution is practically identical 
in form to the nodal BEC wavefunction with two nodes 
inside the well, except for the asymptotic behavior out- 
side the well — see the solid and dashed lines in Fig. 0] 
[2l| . Inside the well the nonlinearity in Eq. ^ has little 
effect on these nodal solutions and the nodal BEC wave- 
function strongly resembles single-particle orbitals. To 
stress the similarity, we will refer to the BEC solutions 
discussed above (the solid lines in Fig. [3]) as the "atomic 



nodal BEC 

NI E = scattering solution 
NI bound state (1 node) 
NI bound state (ground) 




FIG. 4: For a well depth of 77 = 3.5, the noninteracting (NI) 
bound states with zero and one node and the NI zero energy 
scattering solution with two nodes are plotted to compare to 
the interacting solutions with the same well depth in Fig.[3jc). 
In particular, the (green) nodal BEC solution with two nodes 
is similar to the zero energy scattering wavefunction. 



nodal BEC" solutions. 

We show that the wavefunction of the other solutions 
(with a node that moves further away from the well edge 
as the well-depth increases), plotted by the red and blue 
dashed lines in Fig. [3]Jb-c), are related to the nonlinear 
physics of solitons. Note that the curves of the same color 
in Fig. [2Ib-c) have the same peak magnitude and would 
have similar shapes (except for the outermost node) . We 
label the solid line solutions as -i/'a.i [f) and the dashed so- 
lutions as '!/;n.i(r), representing the type of solution by the 
first ("atomic nodal BEC" and "nonhnear nodal BEC") 
and the number of nodes by the second subscript. The 
"nonlinear nodal BEC" solutions, when they exist, are 
well approximated as 



V'n,j(^) ~ ■4^a.,i-i{r)x{r ~ R), 



(4) 



where x(r) = tanh(77^) and R is the position of the 
outermost node. This hyperbolic tangent is the familiar 
nonlinear kink or soliton solution: the ipn.i solution has 
a soliton that is trapped or pinned outside the potential 
well. We illustrate this point graphically in Fig. [5] by 
plotting tpn,i{r)/'4'a,o{'r) for different well depths 77. Note 
that the quotient follows the hyperbolic tangent shape. 
The red dashed (1 node) wavefunction in Fig. [3] can be 
interpreted as the red solid (nodeless) wavefunction with 
a kink superimposed on it. Similarly, the blue dashed 
(2 node) wavefunction can be viewed as the blue solid 
wavefunction with a kink. For deeper wells, this pattern 
repeats, e.g., for deeper wells the green (2 node) wave- 
function develops a kinked partner with 3 nodes. 

Which solution has the highest energy? By fixing the 
density far away from the localized potential we fix the 
chemical potential so that we have to determine the free 
energy, F = H — fj,N, where H is the Hamiltonian and 
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an attractive well smoothly lowers AF from to negative 
values as the well is turned on. At the well depths where 
the nodal solutions first appear, their AF is positive for a 
range of well depths before becoming negative [2^ . The 
nonlinear pinned soliton solutions of n nodes (depicted by 
the red and blue dashed lines in Fig. [6]) have a AF value 
that differs from that of the atomic n — 1 node solutions 
by a nearly constant value — the energy cost of creating 
the kink. 

The solutions with nodes are obviously not the ground 
state of the system, but are they metastable? In the next 
sections, we undertake a stability analysis to investigate 
whether they might be long-lived. 



FIG. 5: Plots of '4^n,i{r)/'4^!!,fi{r) for various well depths; the 
quotient takes on the hyperbolic tangent form discussed below 
Eq. Q, indicating that the solution tpn,! can be interpreted 
as a kink added to the solution i/'a,o- 



N ~ J d^r \tpo{r)\'^ ■ In fact, we calculate the difference 
between the free-energy of the <f>o solution and that of 
the homogeneous system of the same density (without 
external potential), 



AF 



|V<^o| 



1 



+ Uo\<f>o\' + -i\M' + ^)-^i 



(5) 

In Fig. [HI we plot AF vs. the well depth rj for the nodeless 
and the two nodal solutions for a well radius of 6 = 1.0 ^. 
For ease of comparison, the colors/dashings of the free 
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FIG. 6: Free energy difference between the system with the 
potential well present and the homogeneous system without 
well, as given by Eq. [5] for a well radius of & = 1.0^. The 
colors/dashings correspond to those in Fig. O 

energies in Fig. [6] are the same as those of the wavefunc- 
tions in Fig. [31 

As with noninteracting systems, the nodeless solution 
(red line) has the lowest free energy. The introduction of 



III. LINEAR STABILITY ANALYSIS: 
FORMALISM 

In this section and the following we test whether or 
not the (jjo function is stable with respect to small per- 
turbations. We use linear stability analysis, which inves- 
tigates the response of the system to a weak perturba- 
tion caused, for instance, by a modulation of the poten- 
tial around the time-independent external potential Uo, 
U (x, t) = C/o(x) + (x, t). The response of the system is 
described by the time-dependent Gross-Pitaevskii equa- 
tion. 



. g0(x,t) 

' dt 



-^ + |(/.(x,i)|2 + [/(x,0 



,t). (6) 



The time-independent solution to the equation with U ~ 
Uo evolves as (/>(x, i) = exp{—it)(po{x.). The system's 
small amplitude response 



(x, t) « cxp(-it) [(?!)o(x) + 6(j>{x, t)] , 



(7) 



to the potential perturbation U{x,t) ~ [/o(x) + SU{x,t) 
evolves according to the linearized Gross-Pitaevskii equa- 
tion for the 66 evolution. 



I—- = ft bcp + t 
at 

where the h~ operator, 
h- 



' + (f>oScf>*) + (f>oSU, (8) 



^ + C/o+(|0op-l) 



(9) 



represents the Hartree-Hamiltonian of the (jjo system. 

Note that the linearized Gross-Pitaevskii equation, 
Eq. ([8]), couples 6(j) and d(f>* , mixing wavefunction fluc- 
tuations and their complex conjugates. To describe the 
coupling, we separate out the positive from the negative 
frequency components in the Fourier transform. 



(5(/)(x,a;) = 



dte" 



(x,t); 



dtoe 



(,x,wj 



(10) 
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Wc define the positive, 5(j)p and negative, <50„, frequency 
amplitude functions to have non-zero-values for w > 0, 

(50p(x,a;) ~ 6(f){-x.,uj) if w > 0; 

= if w < 0; 

50„(x,a;) = (5(/>*(x, — cj) if w > 0; 

= if < 0: 



so that 



(50(x, u) = S(f)p{x, uj) + S(f)*^{x, —uj). 



(11) 



(12) 



If the time-dependent function (5[/(x, t) is real-valued, 
then the negative and positive frequency amplitudes are 

identical, SUn(x,uj) = SUp(x,uj), where 6Un and SUp are 
defined as in Eq. In the time-domain. 



Aujt 



■(?!)p(x,i) +(5(?i;(x,t) 



(13) 



where, in accordance with the Fourier transformation, 

<5(/)p(„)(x,t) = — / e (5(/)p(„)(x,cj), (14) 

and we remember that (50p(„) (x, tj) = if w < 0. 

Inserting Eq. into the linearized Gross-Pitaevskii 
equation of Eq. wc identify the e~*"* and the e*"* 
components. Taking the negative of the complex conju- 
gate of the e*'^* components, we obtain 



\<f>o\ 



= -4'l*S(t)p - 



h- 



<j>oSUp; 

, - ^oSU^- (15) 



Multiplying by e and integrating over (27r) ^ J duj, 
we obtain the time-domain version, 



d 



C 



u^p 



4>oSUp 



(16) 



where we have introduced the C operator used in Ref. |23l . 

a, 



c 




(17) 



The equations that diagonalize the £ operator, 

£( -'r:_( I - co, ( -'y:_( i , (is) 



Wj(x) 



Uj(x) 

«j(x) 



are the Bogoliubov-de Gennes (BdG) equations [23 . 

If LOj is an eigenvalue of C then lo* is also an eigenvalue. 
If (f>o is real-valued, this statement can be verified by 
taking the complex conjugate of the BdG equations, but 
this property is true for complex- valued (/)o wavcfunctions 



as well. Note that if 



i'j(x) 



is a right-eigenvector of 



£, as shown in Eq. (HH]), then (it*(x), -u*(x)) is a left- 
eigenvector of eigenvalue to* . By integrating 



/ 



d'- (-Mx),-.Mx)) ci^lf^^ 

Uj) ~ {Vj>\Vj)) , 



(19) 



where {uji\uj) — J d'^.Tu*, (x)w*(x), wc find that if u*, ^ 
^j^ {uj'Wj) - {vj'\vj) = 0. 

The above orthogonalization relation suggests the 
{{u\u) — {v\v)) form as a scalar product and we could try 
to normalize the (uj,Vj) solutions by requiring {uj\uj) — 

Mj-+(x) 



1. This, however is not possible: if 



is a right eigenvector with eigenvalue uij , then direct sub- 



stitution shows that 



_(x)\ /z;* + (x) 



eigenvector of eigenvalue 



If the 



is a right 



is normalized by requiring 

then (ui -\ui -) — (vi -\vi -) = — 1. 



j^j^ eigenvector 

- {Vj,+ \Vj,+) = 1 

j,j^_|Lij^_/ — \iyj^_|iyj^_/ — It is then conve- 

nient to distinguish between a "-f" family, — 
{vj\+\vj^+) = 6jij and a " family with {uji ^^\uj^-) — 

We expand the positive and negative frequency com- 
ponents of the wavefunction fluctuations as 



ip(x,i) ^ _ 
i>„(x,t) 



Uj,s(x) 



(20) 



where the s subscript indicates the family, s = -1-1,-1, 
and the (uj,<i, Vj,s) denotes the right eigenvector of C with 
eigenvalue Wj^g- Substitution of Eq. (|20|) into the lin- 
earized Gross-Pitaevskii equation, Eq. (|16p . yields 



k,s' 



dCk,s'{t) ( Mfe^,,(x) 



dt 



"fc,s'(x) 
Vk,s' (x) 



(/)o(x) (5C/p(x, t) 

-0*(x) <5f/„(x,i) 



(21) 



Multiplying from the left by (u 
over the position coordinate x wc obtain 



■J, SI 



3.S, 



and integrating 



9t 



d^a: u* ,,(x)(/)o(x)5[/p(x, t) 



d3xi;*,(x)0:(x)<5t/„(x,t) 



(22) 



Assuming that the perturbation was turned on at a time 
to, the solution that satisfies the corresponding boundary 
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condition, Cj^s{t') = when t' < to, is 



to 



X J <Px' [u,-,(x')0o(x')'5f/p(x',r) 
+z;,-,(x')0:(x')<5C/„(x',r)]. 



(23) 



If LOj^s has a finite imaginary part, one of the c amph- 
tudes grows cxponentiaUy. Even if the imaginary part 
of the (j, s) frequency corresponds to a damped mode, 
ujj^s = Lo — iV, another mode (n, s) exists for which 
LOn^s = s — ^ + «r, so that c„_s(<) grows exponentially. 
In that case, even the smallest of perturbations gives rise 
to an exponentially growing response. In physical and 
simulated BECs, the effects of perturbations are always 
present, caused either by thermal or quantum fluctua- 
tions in physical systems or by round-off errors in nu- 
merical simulations. From these considerations the cri- 
terion for dynamical stability in linear stability analysis 
follows: all eigenvalues of the BdG equations need to be 
real-valued for the system to be dynamically stable. 

For studying the stability of the standing wave Gross- 
Pitaveskii solutions of Section II (in which case the (f)o 
wavefunction, and hence C, is real-valued, and 0^ ~ 
|0oP)j we find it convenient to calculate the sum and 
difference vectors instead of the u and v functions. 



/ 



(24) 



where, from now on, the subscripts j fully characterize 
the mode (including the "family"). By adding and sub- 
tracting Eq. p8)) the BdG equations take the form of a 
diagonal set of coupled eigenvalue-type equations. 



(25) 



where 



1 d d 



2p dp dp 2p' 

I d d iri^ rr r \ 

-2-pYp%+2^ + ''°^P^ 



1; 



(28) 



denote the radial operators for fixed azimuthal quan- 
tum number. 

What is the physical meaning of the /+ and /~- 
fluctuation functions? At finite oscillation frequency, 
is proportional to the density fluctuation and /~ to the 
phase fluctuation. This can be seen by writing the wave- 
function as (/) = -v/ne*", defining the equilibrium density 
and phase values as rio, ao, and introducing their fluctu- 
ations n = rio + 6n, a = + Sa, so that 



6n 



(29) 



We see that we can isolate the density and phase fluctu- 
ations by taking 



5n 



~ 2is/riZ5a. 



(30) 



But from Eq. (|13p . 5(f> ~ 6(f)p + 5(j)^, and putting 5(j)p = 
Uj and 5(j)n — Vj if only the j-mode is excited by the 
perturbation; we find, with Eqs. ([24|) and ((26|) . that 



/+(p)e™«e— * 

^im9 ^ — iuJrnt 



5(j)~5(j}* = /„(p)e^ 
Comparing Eqs. ((30|) and (fSTI) . we obtain 



c.c. 
c.c. 



(31) 



dn — 2y^?i^Re 
1 



5a 



:Im 



/™(p)e™''e 



(32) 



where the h operators, 
2 



[/o + [2cf,l - 1] ± 



are single-particle-Hamiltonian-like operators. For a 
cylindrically symmetric potential and standing wave pat- 
tern, the modes will be eigenvectors of the angular mo- 
mentum operator. We perform the separation of vari- 
ables into radial and angular functions. 



/f(x) 



ftApy 



me 



(26) 



where p is the dimensionless radial coordinate measured 
in units of the healing length ^, 6 is the azimuthal angle, 
and |m| ~ 0, 1, 2, the azimuthal quantum number. 
Henceforth we will take the radial quantum number i> as 
implicit. Using Eqs. (|24p and (|26p . the BdG equations 
take on the form 



mJ ni 
+ 



^rafrr 



(27) 



so that and are indeed proportional to the density 
and phase fluctuations, respectively. 

The coupled BdG equations, Eqs. (P7|) . can be cast in 
an alternate form by noting that the repeated application 
of LJrn yields decoupled eigenvalue equations. 



m J rn 
,2 



h~ h+ f+- 
iJ ri 



(33) 



so that /+ is the eigenvector of h^^h^, whereas /~ is the 
eigenvector of h^h:^. While this eigenvalue problem is 
distinct from the Hamiltonian eigenvalue problem {h^h^ 
has a fourth derivative in space, for instance), there is one 
important similarity: since the product operators h^h^ 
and h~^h^ are both self-adjoint, their eigenvalues are 
real- valued. As a consequence, ujm is either entirely real 
(if > 0) or entirely imaginary (if < 0). The former 
case, real uJrm corresponds to oscillations about a possi- 
bly stable equilibrium solution (j)o for mode to; whereas 
an imaginary ujm indicates that the mode to drives the 
system away from an unstable equilibrium solution 0o- 
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IV. LINEAR STABILITY ANALYSIS: 2D BEC 
STANDING WAVE PATTERNS 

We subject the standing-wave 0o wavefunctions of Sec- 
tion II to the stabihty analysis of Section III. We use 
the form of the BdG equations, Eq. ((28|) . and solve 
Eqs. ([33|) for and f^. In the following two subsec- 
tions we describe typical examples to uncover the gen- 
eral stability picture. The first case deals with a shallow 
well that supports nodal BECs with a single node and 
the second case deals with deeper wells and nodal BECs 
with many nodes. We highlight the following key results: 

(i) the nonlinear nodal BECs are always unstable, and 

(ii) the atomic nodal BECs can transition from stable to 
unstable and from unstable to stable as the well depth 
increases. 



1. Single node BEC states in a well of coherence length 
radius 

We consider the stability of a BEC with one node in- 
duced by a potential well of radius b ~ 1.0^, as pictured 
in Fig. O While the azimuthal quantum number m has 
to be an integer to ensure the single-valuedness of /,„ in 
Eq. we can formally solve Eqs. (P7|) for any value 

of m . Varying m continuously, we plot Im[cj„i] for the 
nonlinear (a) and atomic (b) nodal BEC solutions with 
one node for various well depths. 

When the atomic and nonlinear nodal solutions first 
appear with increasing 77, they are identical. For a well- 
depth of 77 = 1.25, we show the two wavefunctions with 
one node in Fig. ^a.). The blue line shows the atomic 
nodal BEC wavefunction and the red dashed line rep- 
resents the nonlinear nodal BEC wavefunction. Note 
that in spite of their similarity, the Im[wm] spectra differ 
markedly, as seen in Fig. [71 The nonlinear nodal BEC 
has a finite value of Im[ci;„i=o]i while the atomic nodal 
BEC has Im[ci;m=o] = 0, a general features that distin- 
guishes atomic from nonlinear solutions. In Fig. [71 we 
see that a small variation of ?/ from rj ~ 1.2445 (within 
0.0001 of the critical well depth for nodal solutions) to 
77 = 1.25, drives the left-edge of the Im[wm] lobe in oppo- 
site directions: Im[u;m=o] = for the critical well depth 
where the nodal solutions first exist; and uim=o becomes 
entirely imaginary for the nonlinear nodal BEC solutions 
and entirely real for the atomic nodal BEC solutions. 
The nonlinear nodal BECs have at least one imaginary- 
valued frequency mode, the radially symmetric m — 
mode, so that the nonlinear nodal BECs are always un- 
stable. The instability of the nonlinear nodal BECs has 
features in common with the snake instability of a kink- 
wise nodal line or plane in an otherwise homogeneous 2D 
or 3D condensate but an in-depth discussion would 
take us further afield. 

We now discuss the atomic nodal solutions. Fig. ^h) 
shows Im[ti;m] for the one-node atomic nodal BEC state. 
As mentioned above, Im[ti;m=o] = and the lowest possi- 



' I ' ' ' ' I ' ' ' ' I 
(a) nonlinear, 1 node 




G-^. ri= 1.2445 
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FIG. 7; The linear stability spectrum Im[ci;m] is shown for 
several well depths 77 for (a) the nonlinear nodal BEC and (b) 
the atomic nodal BEC, both with one node for well radius b — 
1.0^. The noninteger value jt^'s and the connecting lines have 
been included to guide the eye, although only integer-valued 
azimuthal quantum number m 's are physical. The nonlinear 
nodal BEC (a) always has at least one integer m with an 
imaginary Um and thus is always unstable. The atomic nodal 
BEC (b) has no imaginary-valued uim for integer m when 
77 = 1.8, 2.2, 2.6, and 3.0, and is thus stable for these cases. 



ble angular momentum of an exponentially growing mode 
is 777 > 1. For the lowest three well depth values shown 
in Fig. [7{b), the m = 1 mode has an imaginary- valued 
frequency and the atomic nodal BEC is unstable. As the 
well increases, from 77 = 1.8 to 77 = 2.2, for instance, we 
don't find any imaginary values for integer m and 
the atomic nodal BEC is stable. Between 77 = 2.2 and 
77 ~ 2.6, the Im[ti;„i] lobe shifts to larger m values and 
passes through the m = 2 line, indicating the instabil- 
ity of the mode with azimuthal quantum number 2. In 
the 77 = 2.6 to rj = 3.0 interval, the Im[tj,„] lobe has 
shifted to values between m = 2 and ttt, = 3, so that this 
regime of well-depths is stable again. Therefore, as the 
potential deepens from 77 = 1.25 to 3.0, the single node 
atomic nodal BEC transitions from instability with expo- 
nentially growing m = \ mode to stability, then returns 
to an unstable status, this time caused by the exponential 
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growth of the m = 2 mode, then enters another island of 
stabiHty. We refer to this trend as reentrant stability to 
describe the entering and exiting into and out of islands 
of stability. 



2. Multi-node BEC states in broad, deep wells 

We now discuss nodal BECs in a wider, deeper well 
with radius b = 4.0^ and depth 77 = 11.5. This well has 
nodal BEC solutions with up to six nodes; Fig.[5Ja) shows 
the five node atomic and six node atomic and nonlinear 
nodal BECs. Fig. El^b) shows the corresponding Im[cL'm] 
spectra. As in the previous subsection, we find that the 
nonlinear solution has an imaginary-valued m = mode, 
while the two atomic solutions do not, following the gen- 
eral m = behavior mentioned above. More striking, 
however, are the spectra's multi-peaked structures which 
we now explain. 

To understand the Im[Lj„i] structure we revisit the BdG 



T" 



equations (|27|) . h^f,"^ 



The frequency tOm is 



either entirely real or entirely imaginary. Therefore, if we 
treat m as a continuous variable and we assume lu„i to be 
a continuous function of m then cUm has to pass through 
zero to switch from real to imaginary values. When uj„i 
vanishes, either h^f^ = or h:^f^ ~ 0, corresponding 
to a or /~ eigenvector of the operator with zero 
eigenvalue. 

Consider the operators of ((28)) which we write as 
h-m = —j^-§^P-§^ + Vm introducing the potentials, 



V^^-^ + Uo{p) 



0o(2±l)-l 



(34) 



that are effective potentials of the single-particlc- 
Hamiltonian-like operators /i^. We show examples of 
for m = in Fig. [S^c) for the two atomic nodal 
BECs of Fig. [HJa). The short-scale, oscillatory structure 
of the potentials stems from the 0o contributions, 
but for high-node BEC states in sufficiently deep wells 
(for which 0^ ^ Uq) the overall shape is dominated 
by the external Uq potential. For example. Fig. [SKc) 
shows for the state with five (six) nodes and the 
potentials themselves support five (six) negative energy 
levels. As m increases, the angular momentum barrier, 
m^/2p^, grows and raises the energy-levels to positive 
values [26[. At a value mi, a V^_^ level crosses the zero 
energy dividing line and h^_^f^^_^ = so that the 
BdG equation becomes an eigenvalue equation with zero 
eigenvalue. This mi value denotes the left-hand side of 
the first peak or lobe in the Im[a;m] spectrum. As m 
increases further to TO2 at which the corresponding V^^^ 
level crosses over, h^^f^^ — 0, the BdG equation re- 
duces to another eigenvalue equation with zero eigenvalue 
and the m2 value denotes the right-hand edge of the lobe 
that started at mi. The pattern repeats as many times 
as there are levels in the V^^q potentials, generating the 
other peaks in the Im[Li;m] spectrum. 



— atomic, 5 nodes 

— nonlinear, 6 nodes 

— atomic, 6 nodes 



(a) 





CM atomic, 5 nodes 
CM nonlinear, 6 nodes 
a-e atomic, 6 nodes 



2 4 6 8 10 12 
azimuthal number, m 
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— atomic, 5 nodes 

— atomic, 6 nodes 



(c) 



FIG. 8; (a) Nodal BEC wavefunctions five and six nodes are 
shown for a well radius b = 4.0 ^ and well depth 77 = 11.5. (b) 
The linear stability spectrum Im[ix)m] for the three solutions 
of (a); as in Fig. [T] only integer m's are physical and the 
noninteger m's and lines are shown to clarify the trends, (c) 
The potentials V^^q of Eq. p4|l are useful for interpreting 
the linear stability spectra, as discussed in the text. 



The above construction shows that the number of 
peaks in the Im[a;m] spectrum is equal to the number 
of negative energy levels of the V^^q potentials. More- 
over, the widths of the Im[a;m] peaks, Am, are related to 
the energy splitting of corresponding states in the V^—q 
and in the V^^q potentials. We illustrate this point with 
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Fig. m^c): as the BEC density (j)^ in the potential well 
is smaller in the six node than in the five node BEC, 
the I^j^^o Kn=o potentials resemble each other more 
in the six node atomic BEC state than in the five node 
atomic BEC state. As a consequence, the difference in 
energy between corresponding pairs of V^^q and V^^^ 
eigenstates is smaller in the six node than in the five node 
case. Therefore, the m value at which the and the V~ 
state reaches zero energy differ by a smaller amount Am 
in the six node than in the five node situation. Gener- 
ally, as the atomic nodal BEC with the maximal number 
of nodes has the smallest BEC density 0^, this time- 
independent BEC solution has the smallest peak widths 
Am. 

The above analysis reveals that common methods for 
bound state physics can determine relevant instability 
properties of standing wave BECs. For instance, we con- 
sider the upper imaginary azimuthal quantum number 
TOc above which all mode eigenfrequencies ojm take on 
purely real values. We estimate rric for the atomic nodal 
BEC state of highest node number contained in a po- 
tential well of sufficient depth to ensure that the BEC 
density satisfies |0oP ^ Uo inside the well (for the node- 
less case, this is obviously not true). As we ramp up the 
m quantum number, the last bound state of the po- 
tential becomes unbound at m = rric. The corresponding 
eigenstate (resonant) wavefunction, which we shall 
refer to as (and we assume Xo to be normalized) is 
a nodeless function roughly centered on the position of 
the potential minimum of V^, i.e., on p = h. At the zero 
energy crossing point, 

Y72 2 

(Xo I - — + f^o + ^ + [4>l - 1] IXo ) = (35) 

Neglecting the kinetic energy term and the [0^ — l] con- 
tributions and estimating the angular momentum barrier 
contribution by m^/25^, we find m^ ~ b^/2Uo, or with 
our notation for the square well potential 

rric^Y' (36) 

which overestimates rric because we neglected the kinetic 
energy and the (1)1 terms. Nonetheless, we expect this ap- 
proximation to be reasonable for deep wells. For exam- 
ple, in Fig. [HI the observed nic of ^ 15 can be compared 
to the TTic estimate of 18.1 from Eq. ([55]) . 

All three nodal BEC-states shown in Fig. [8] are un- 
stable. The five node atomic nodal BEC has imaginary 
modes for m = 2,3, 5, 8, 11, and 14. In Fig. [S] we show 
the density fluctuation pattern of the m = 3 mode, pro- 
portional to /+(p)cos(to0) in accordance with Eq. ([32)1 . 
If the m = 3 mode were the only unstable mode or if it 
were the dominant unstable mode (growing at the largest 
rate = Im(a;in)), the density variation would take on 
the form pictured in Fig.[8]in the initial stage {t < r^) of 
the break-up of the nodal BEC state. The intricate oscil- 
latory pattern of Fig. [5] highlights the role of interference 
in the wave dynamics of the instability. Only when the 



interference fringes "match up" properly to initiate the 
break-up does the instability set in. In addition to the re- 
striction of m values to integers, which is an expression of 
the angular boundary condition S(f){p, 9) = 5(/)(p, 9 -I- 27r), 
the radial oscillations have to line up to create the con- 
ditions at which the instability can set in. We interpret 
the finite Am width of the Im[ti>,„] lobes as the range of 
m values for which the radial interference fringes match 
up sufficiently well to trigger the break-up of the nodal 
BEC structure. 



6n [arb.] 




m 

FIG. 9: The density fluctuation for the m = 3 mode, 5n 
from Eq. (|32|l . is shown for the five node atomic nodal BEC 
state shown in Fig. EJa). This is an unstable mode with an 
imaginary ujm and the pattern seen here will contribute to the 
break-up of this nodal BEC. 

The six node atomic nodal BEC whose Im[a-'m] spec- 
trum is shown in Fig. [S] has unstable modes for m — Q 
and m = 15. However, the widths of its Im[a;ni] lobes 
satisfy Am < 1. This means that as the depth of the po- 
tential well is varied, the imaginary-valued ujm intervals 
can shift off of all integer m's leading to a reentrant sta- 
bility behavior similar to that discussed in the previous 
subsection. 

Can we make a general statement about the stabil- 
ity of stationary multi-node BEC wavefunctions? The 
multi-lobed structure of the Im[<x'm] spectrum can actu- 
ally promote stability: if the unstable m intervals. Am, 
are sufficiently narrow, Am ^ 1, the Im[wiii] lobes may 
not overlap any integer m values. The progression from 
the five node to the six node Im[wni] spectrum shown 
in Fig. [S^b) indicates an interesting trend: while the in- 
crease in nodes increases the number of lobes, their l^m 
widths decrease markedly. Furthermore, the Ajti widths 
decrease for wells of smaller radius, as we now show. 

We can derive a rough estimate for the peak widths 
Am along the lines of the above argument that asso- 
ciates the vanishing of LOm with the zero-energy crossing 
of bound state levels of . We assume that the BEC 
density inside the well remains sufficiently small to de- 
termine the splitting between the and eigen en- 



11 



ergies by perturbation theory. The necessary condition 
reduces to \Uo — l\ in the relevant well region — the 

region in which the the Vj^ eigenstates have significant 
density. We determine the and V~ eigenstate en- 
ergies as perturbations on the eigenstate of the average 
potential ^ [V+ + V,-] /2, 



F™(p) = C/o(p) + 202(p)_i + 



2p2 



(37) 



We denote the eigenstate of the Vm potential with i/ 
nodes by Xm,i^ and we denote its energy by E^{m). In 
first order perturbation theory, the energy of the cor- 
responding eigenstates of the and V,~ potentials, 
(m) and E~ (m) , are then equal to 

E^ (m) ^ E^ (m) ± A£;^ (m) , (38) 

where lS.Ei,{m) is the first-order 0^ energy splitting, 

A£;,(m) = (xm,.|0olXm,.> • (39) 

The slope of the m variation of E,j{m)^ dEn/dm, can 
be calculated with the HcUmann-Feynman theorem (the 
derivative of the hamiltonian expectation value with re- 
spect to a hamiltonian parameter is the expectation value 
of the derivative of the hamiltonian operator) : 

^ - ly ) 



P 



(40) 



Then, the difference Atoi/ in m values at which the E^ 
and the values cross the zero-energy line can be cal- 
culated from the matrix elements for m = mo,u at which 
the ly-level of the average Vm potential crosses the zero 
energy line, 

2AE^{rno.^) 



dEy/dm 

(Xmo,J0olXm„ 



(41) 



Assuming that the Xm^ „ function is centered on the po- 
tential edge p = b, we can estimate the denominator as 
TOo,i//6^- Denoting (xrno,J0olXrno,J by , we find 



A? 



(42) 



Since the Xm^ v function is centered on the edge of the 
Uo potential where the (po function is reaching towards its 
asymptotic value, we try an even cruder approximation 
^ 1, resulting in 



2&2 

Ami/ ~ — — 



(43) 



which overestimates the Am widths. Nevertheless, we 
believe that the general trends predicted by Eq. (|43p are 
correct. 

From Eq. we expect that the peak widths Am of 
the Im[cL'„i] spectrum narrows (i) as the lobe is situated 
at larger m values, and (ii) as the radius of the potential 
well decreases. The former effect can be observed in the 
spectrum of Fig. [5]^b) . The latter effect implies that nar- 
row wells with a subcoherence-length radius, 6 <C , have 
peaks of widths Am ^ 1 , implying a decreased likelihood 
of lobe-overlap with an integer m value (hence, increased 
likelihood of stability). 



V. DYNAMICAL STABILITY FROM 
SINGLE-PARTICLE HAMILTONIAN-LIKE 
EIGENSPECTRA: RECIPE 

The treatment of Section IV makes explicit use of the 
cylindrical symmetry by introducing the device of a con- 
tinuously varying azimuthal quantum number m. In this 
section, we propose a different, more general procedure 
to test the dynamical stability of localized, potential- 
induced standing wave BEC patterns. The recipe we 
describe below is based on the same concept of continu- 
ously varying frequencies Uj that have to vanish as they 
cross over from real to imaginary values. At the crossing 
points, the BdG equations reduce to a single-particle- 
Hamiltonian-like eigenproblem for zero eigenvalue. As 
in Section IV, the general recipe allows one to infer the 
stability of a standing wave pattern from the eigenspec- 
tra of single-particle-like operators alone, circumventing 
the need of conducting the full linear stability analysis. 
Unlike the technique of Section IV, the recipe is indepen- 
dent of the symmetry of the potential and of the standing 
wave density \4>o\'^ and can be applied in three as well as 
two dimensions. 

Recipe: We consider the eigenspectrum of the single- 
particle-Hamiltonian-like operators 



V2 



+ f/o(x) + 2 I 



1 ± 



(44) 



"similar" eigenstates j of the 
operators with eigenvalues E'^ and EJ respec- 

indicat- 



We identify pairs of 
and h~ 

tively and with the same quantum number j, 
ing eigenstates that would go over into each other in 
the limit that (/)o ^ in Eq. (|44|) . We suggest that 
the standing wave solution (j)o of the time-independent 
Gross-Pitaevskii equation is unstable if there is one or 
more quantum numbers j for which EJ < and E'^ > 
and that it is stable otherwise. 

Justification : The justification of the recipe is based 
on the structure of the BdG equations and the continu- 
ity of its eigenfrequencies ujj with respect to a continu- 
ously varying parameter e that appears in the Hamilto- 
nian. Such parameter could be introduced in different 
ways (multiplying the chemical potential or the interac- 
tion strength, for instance). Here we choose to vary the 
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depth of the external potential 

C/e(x) =eC/o(x), 



(45) 



and we let e range from to 1 , assuming that a standing 
wave solution has appeared by e = ec, with ec < 1- The 
standing wave function 0o,e solves the time-independent 
Gross-Pitaevskii equation 



+ C/e(x) + |(/.o,e(x)p-l 



0o,.(x) = (46) 



with e > ec- The BdG equations for fixed e and real- 
valued 0o function, 



(47) 



with 



ht = -— + [/e(x) -f 20o,(x)2 - 1 ± 02^^(x) , (48) 

has solutions that are pairs of functions {f^^,f^^)- As 

an eigenvalue of the Hermitian h'^h^ operator, ujj{e) is 
real- valued so that loj (e) is either entirely real or entirely 
imaginary. To cross from real to imaginary values, ujj{e) 
has to pass through zero, at which point the BdG equa- 
tions reduce to a Hamiltonian eigenvalue problem of van- 
ishing eigenvalue. If xj is the zero-eigenvalue, the j-type 

eigenvector of the h~ operator for e = e~ , 



dividing line by e = 1 whereas > 1 so that E^' > 

and the hj' eigenvalue has not been brought below the 
zero-energy line yet, then ujj has moved through only 
once, switching its value from real to imaginary value. 
If both E~ < and < 0, then ijjj{e) moved twice 
through as e varied from ec to 1 switching from real to 
imaginary and back to real. The upshot is that if both 
eigenvalues E^ find themselves on the same side of zero, 
then u!j is real, if they are on opposite sides of the zero 
energy dividing line, tOj is imaginary. 

The application of the above method to the 2D- 
cylindrically symmetric situation described in Sections 
II and IV leads to the investigation of E^ for integer, 
fixed values of m = 0, ±1, ±2, from which stability or 
instability can be determined. In contrast, the treatment 
of Section IV used the quantum number m itself as the 
continuously varying parameter and calculated the rate 
of the instability Im[w,„]~^ as a function of m. A single 
Im[a;„i] curve then reveals excitations of different j-types: 
the different Im[a;m]-lobes shown in Fig. ^ correspond 
to /-functions of different j-types, here modes of different 
vibrational quantum numbers u (with different numbers 
of nodes). While convenient for studying the 2D cylin- 
drically symmetric standing wave pattern, the method of 
a continuously varying quantum number cannot readily 
be generalized to describe non-symmetric 0o-patterns or 
three-dimensional structures; but the concept outlined 
here of varying a parameter in the effective Hamiltonian 
can be. 



(49) 



VI. CONCLUSIONS 



then {f^ = 0, fj ~ Xj ) solves the BdG equations for 
that specific e value. Conversely, if x'j' is the zero- 
eigenvalue, j-type eigenvector of the hf operator for 



K.xt-0 



(50) 



then (/+ 



0) solves the BdG equations for 



e — e^. In general, hfxj.e = ^fi^)xf,e ^^"^ the zero 
eigen energy crossings happen at specific e- values, e = 
for which Ef{ef) = 0. As (0,xj) and (x^,0) solve the 
BdG equations for specific values of the parameter e, the 
functions must be of the same type (same number 

of nodes, same quantum numbers) as the xf wavefunc- 
tions. As uij{e) passes through zero at e = e^, it either 
changes sign or switches between real and imaginary val- 
ues. We argue for the latter: as an eigenvalue of the 
hfh~ operator, ^j(e) and its derivative with respect to 
e are continuous. As a consequence, w|(e) switches sign 
at e = , corresponding to tUj (e) switching between real 
and imaginary values. If < 1 so that the E~{e) eigen- 
value has already been lowered through the zero-energy 



We have studied the existence and dynamical stabil- 
ity of stationary BEG states with nodes in the BEG 
wavefunction. In particular, we investigated large 2D 
BEG systems in the presence of an attractive potential 
well with a coherence length-sized radius. The potential 
is cylindrically symmetric and of the square well type. 
From the center of the potential well, the BEG tends to 
a uniform density. As the depth of the attractive well in- 
creases, BEG solutions with radial nodes in the wavefunc- 
tion appear in pairs, each pair exhibiting one more node 
than the last pair. While the calculations are carried 
out for this very specific case, we discuss which trends 
are independent of the symmetry and can be applied to 
three-dimensional standing wave BEG patterns. 

We have classified the nodal BEG solutions to the 
time-independent Gross-Pitaevskii-equation to be of two 
types. One type we refer to as the atomic nodal BEG 
state since the wavefunction resembles that of single- 
particle-like bound states in the attractive potential well. 
The other type we refer to as the nonlinear BEG state 
since the n node BEG wavefunction resembles the n — 1 
node atomic BEG wavefunction with a soliton attached 
outside the potential well. 

Our stability analysis shows that the nonlinear nodal 
BEG state is always unstable and that atomic nodal 
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BEC states can transition from stable to unstable to sta- 
ble status as the attractive well depth is adiabatically 
increased — a behavior which we refer to as reentrant sta- 
bility. In general, the stable regimes grow as the radius 
of the attractive potential is decreased below the BEC 
healing length scale. 

We describe a more general analysis of the dynamical 
instability that is based on the cigcnspcctra of two ef- 
fective single-particle-Hamiltonian-like operators, which 
gives new insights, can be applied to low symmetry stand- 
ing wave patterns patterns, and is valid in three dimen- 
sions. 

Given the dynamical stability in specific parameter 
regimes, atomic nodal BEC states could find applica- 
tions, for instance, in a BEC ring as a rotation sensor. 
Nearby instabilities may be exploited as a threshold mea- 
surement tool. 

This work serves as an introduction to the existence 
and dynamical stability of nodal BEC states. Many chal- 
lenges, on the theory and on the experimental front, re- 
main. We have shown that nodal BEC states can be dy- 
namically stable, but the manufacturing of these states 
may require creativity. Perhaps the phase imprinting 
techniques that have generated vortices and dark soli- 



tons [27|, |28| can create the required initial conditions. 
Conversely, the non-adiabatic dynamics of BECs in steep 
localized potential wells may exhibit signatures of the in- 
stabilities described above: if the BEC dynamics brings 
the system close to a nodal wavcfunction, it can break 
up if the potential depth is such that the nodal state is 
unstable. 

We believe that standing wave BECs with nodal sur- 
faces separated by macroscopic scale distances suggest 
new prospects and we hope that this work will stimulate 
interest in standing wave BECs. 
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